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Biopolymer Networks play an important role in coordinating and regulating collective cellular 
dynamics via a number of signaling pathways. Here, we investigate the mechanical response of a 
model biopolymer network due to the active contraction of embedded cells. Specifically, a graph 
(bond-node) model derived from confocal microscopy data is used to represent the network mi¬ 
crostructure, and cell contraction is modeled by applying correlated displacements at specific nodes, 
representing the focal adhesion sites. A force-based stochastic relaxation method is employed to 
obtain force-balanced network under cell contraction. We find that the majority of the forces are 
carried by a small number of heterogeneous force chains emitted from the contracting cells. The 
force chains consist of fiber segments that either possess a high degree of alignment before cell con¬ 
traction or are aligned due to the reorientation induced by cell contraction. Large fluctuations of 
the forces along different force chains are observed. Importantly, the decay of the forces along the 
force chains is significantly slower than the decay of radially averaged forces in the system. These 
results suggest that the fibreous nature of biopolymer network structure can support long-range 
force transmission and thus, long-range mechanical signaling between cells. 

PACS numbers: 87.15.rp, 87.85.jc 


1. INTRODUCTION 


The extracellular matrix (ECM) is an interconnected 
network of biopolymers that provides structural support 
for cells and allows the diffusion of biochemicals within 
tissues. The most abundant component of ECM is type I 
collagen, a fibrous protein responsible for giving the ECM 
its material stiffness jT] . Cells attach and move through 
the ECM using protein complexes that link the ECM to 
the force-generating cell cytoskeleton Q. However, these 
cell-ECM adhesions also act as sensors, sending infor¬ 
mation to the cell about the structure and mechanical 
properties of the surrounding matrix 4} and helping 
to regulate cell behavior such as motility, morphology, 
and differentiation BH- The stiffness and the relative 
alignment of fibers in the network are particularly im¬ 
portant to cell function. For example, dense and rigid 
collagen gel can promote growth and progression of can¬ 
cer cells and tumors [§]. Other important examples 
are durotaxis in which cells tend move in the direction 
of increasing matrix stiffness 0, and contact guidance 
in which cells tend to align and move in the direction of 
fiber alignment 00 . 

Cell-ECM interaction is a dynamic process in which 
the cell actively remodels the network [3, 0] and these 
effects can propogate over long distances, even affecting 
the bulk propeties of the network 0, 0 ■ Specifically, 
tension exerted by the cells can align the fibers in the net¬ 
work leading to long range force transmission 000. 
Fiber mediated stresses can trigger mechano-sensitive 
pathways of distant cells affecting behaviors such as force 
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generation 00 and cell-ECM adhesions 0 ] and lead¬ 
ing to diverse collective behaviors [2lj| . This coupling of 
cells provides a means for mechanical communication and 
plays an important role in regulating and coordinating 
collective cellular dynamics in a wide range of biophysi¬ 
cal processes, such as morphogenesis, tissue regeneration, 
and immune response, as well as diseases such as muscu¬ 
lar dystrophy, fibrosis, and cancer pH. 0. 0 [2bl | . 


Due to their effect on cell behavior and communication, 
a significant amount of work has been carried out to char¬ 
acterize the structural and physical properties of biopoly¬ 
mer networks. Traditional morphological descriptors for 
such networks include the distribution of fiber length 0] , 
porosit yMJl, pore-size distribution 0-0 and tur¬ 
bidity |33[ H , which are mainly bulk averaged proper¬ 
ties. Recently, local topological and geometrical statistics 
such as distribution of number of fibers at a cross-link 
node (i.e., valency number) and relative fiber orienta¬ 
tions (i.e., direction cosine) are employed to successfully 
reconstruct type I collagen (COL-I) network computa¬ 
tionally 0 ]. Higher order spatial fluctuations has also 
been utilized to characterize the evolution of COL-I net¬ 
work during gelation process 0]. In addition, the trans¬ 
port properties (e.g., macromolecule diffusivity) 0,Ez- 
1461 1 an d mechanical properties (e.g., elastic modulgbulk 
rheology, stress distribution, etc.) fTdl. 0.0. l35l |47I453| 
of biopolymer network, which are respectively crucial to 
the chemical and mechanical signalling between the cells, 
strongly depend on the network microstructure. 


In most studies, the biopolymer network is treated as a 
material system with no biological cells embedded. How¬ 
ever, it is known that cellular interaction with COL-I 
can induce dynamic remodeling of the network. These 
cellular scale dynamics propagate to much longer range 
and affect the bulk properties of the network 0], as well 
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as the behavior of other cells through mechano-sensitive 
pathways and machineries 0, [24], Hg . In addition, fiber 
mediated stress-regulated force generation 0, [3 and 
stress-reinforced cell-ECM adhesion 0] provide internal 
feedback mechanisms that support the emergence of di¬ 
verse collective behaviors jUJ. As a result, understand¬ 
ing the homeostasis of cellularized biopolymer network is 
an essential step towards the understanding of complex, 
self-organized multicellular dynamics [23j] . 

Recently, continuum models of cellularized collagen 
network have been developed to investigate the fiber- 
mediated mechanical coupling between the cells 0,111]. 
Specifically, high-resolution 2D confocal images of cel¬ 
lularized ECM are thresholded to generate a three-phase 
heterogeneous system including a cell phase, a fiber phase 
and an interstitial fluid phase [0 . Finite element anal¬ 
ysis is then employed to obtain stress distribution in the 
ECM. It has been shown that explicitly considering the 
flbreous nature of ECM is necessary to capture the force 
transmission between cell pairs |l6( and the invasion of 
cancer cells in ECM with well defined geometry [H3 |. 
Such behaviors can not be correctly reproduced using ei¬ 
ther linear or nonlinear homogeneous ECM models. De¬ 
spite all these insights, the continuum models do not ex¬ 
plicitly incorporate cell induced fiber re-orientation and 
other structural remodeling mechanisms for cellularized 
ECM, and thus might not accurately reveal the pathways 
for force transmission in the system [0 . 

In this paper, we systematically investigated the me¬ 
chanical behavior of cellularized ECM, in particular, 
fiber-mediated transmission of forces generated by active 
contraction of embedded cells. The biopolymer network 
is represent by a graph (i.e., bond-node) model derived 
from confocal microscopy data 0]. The cell contraction 
is modeled by applying correlated displacements at spe¬ 
cific nodes (representing the focal adhesion sites), e.g., 
displacing the nodes towards a common center. A force- 
based stochastic relaxation method is employed to obtain 
force-balanced network under cell contraction. We find 
that the majority of the forces are carried by a small num¬ 
ber of heterogeneous force chains emitted from the con¬ 
tracting cells. The force chains consist of fiber segments 
that either possess a high degree of alignment before cell 
contraction or those that are aligned due to the reorien¬ 
tation induced by the cell contraction. Fluctuations of 
the forces along different force chains are observed. Im¬ 
portantly, the decay of the forces along the force chains 
are significantly slower than the decay of radially aver¬ 
aged forces. These results suggest that network structure 
could support long-range force transmission and thus, 
long-range mechanical signalling between cells. 

The rest of the paper is organized as follows: In Sec.II, 
we describe our model of biopolymer network and em¬ 
bedded cells. In Sec.Ill, we provide and validate the 
force-based relaxation method that minimizes the strain 
energy of the network to obtain force-balanced network 
configuration. In Sec.IV, we discuss the properties of 
stressed biopolymer network due to the contraction of a 


single cell. In Sec.V, we present results of the mechanical 
response of biopolymer network due to the simultaneous 
contraction of a pair of cells. In Sec.VI, we provide con¬ 
cluding remarks. 


2. NETWORK AND CELL MODELS 




FIG. 1: (Color online). Microstructure of the model biopoly¬ 
mer network studied in this paper (a) and the associated 
pore-size distribution function P(S) (b). The linear size of 
simulation box shown in (a) is 30 pm. 


The biopolymer network of interest is represented as 
a graph, i.e., a collection of bonds and nodes [35]. In 
this model, each bond represents a fiber and each node 
represents a cross link, see Fig. 1(a) [5] The graph rep¬ 
resentation is obtained by thresholding and skeletonizing 
a 3D stack of confocal microscopy images of type I col¬ 
lagen at density of 4.0 mg/ml. The average fiber length 
d° = 1.28pm. A cubic simulation box with a linear size 
of 30pm is used, which contains 5,000 nodes and 8,000 
bonds. The pore-size distribution function of the network 
is shown in Fig. 1(b) [2] 0], Also shown is the Ogston 
approximation for the pore size distribution, i.e., 

P(S) = 2 ( fe(£ + «) e —0 2 (s+ a ) 2 /a ^ (i) 

a z 

where <j >2 = 1 — </>i is the volume fraction of the fibers 
and a is the fiber radius. We note that m is derived 
for networks with very long and stiff fibers, and thus, 
overestimates the occurrence probability of intermediate 
sized pores 0]. 

In order to study the mechanical behavior of the net¬ 
work, we need to specify the properties of individual 
fibers. In this paper, we assume that the fibers pos¬ 
sess short persistent lengths, i.e., the bending modulus 
of the fibers is significantly smaller than the elongation 
modulus. In addition, we consider that cell contraction 
can only generate small forces and thus, the system is in 
the linear elastic regime 0] and the effects of interstitial 
fluid, which quickly dissipates the kinetic energy gener¬ 
ated due to cell contraction, are not explicitly consid¬ 
ered. The mechanically equilibrated network possesses 
the minimal elastic (strain/stress) energy. The elonga¬ 
tion modulus of the fibers is set to be EA = 8 x 10” 7 TV 
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[35j , where E is the Young’s modulus of collagen and A is 
the cross-section area of the fiber. The elastic energy of 
a fiber Ijj defined by nodes i and j) is a quadric function 
of fiber elongation, i.e., 

FA 

e ij = T^cT ' — ^ij) = kj ' (dij — d%j) j (2) 

where the spring constant of fiber kjj = EA/d is in¬ 
versely proportional to its original length . For com¬ 
putational convenience, the coordinates of the nodes in¬ 
side the 30 pm x 30 pm x 30 pm box are rescaled to fit in 
a 1 x 1 x 1 simulation box. In the subsequent calcula¬ 
tions, the unit of length is thus chosen to be the length 
of the simulation box. We note that the results can be 
easily rescaled to the actual units in order to compare 
with experimental measurements. The total elastic en¬ 
ergy associated with the fiber network is given by 

Eo= J2 ( 3 ) 

<i 3> 

where the summation is over all pairs of connected nodes 
(i.e., all fibers lij). 



FIG. 2: (Color online). A schematic illustration of our cell 
contract model. We model an embedded cell as a sphere with 
radius R c (red circle) centered at a prescribed location in 
the network (dark yellow lines). The nodes of the network 
within a certain distance SR to the cell surface (i.e., R c + SR 
to the cell center) are considered as “coarse-grained” focal 
adhesion sites (blue dots). Cell contraction is then modeled 
by displacing the focal adhesion nodes towards the cell center 
with a prescribed magnitude. 

We model an embedded cell as a spherical region with 
radius R c centered at a prescribed location in the net¬ 
work. Although actual cells generally possess much more 
complex morphology, such a simple shape is sufficient for 
our current study, which focuses on heterogeneous force 
chains generated due to cell contraction. The nodes of 
the network within a certain distance SR to the cell sur¬ 
face (i.e., R c + SR. to the cell center) are considered as 
“coarse-grained” focal adhesion sites, through which the 
cell is mechanically coupled with the network [56j|. Cell 


contraction is then modeled by displacing the focal ad¬ 
hesion nodes towards the cell center with a prescribed 
magnitude. In particular, with the boundary nodes kept 
fixed (a node is a boundary node if its distance to the 
box boundary is smaller than 5% of the box length), the 
set of adhesion nodes are displaced towards to cell cen¬ 
ter and then kept fixed. This local perturbation leads 
to a stressed network and the resulting force-balanced 
configuration is obtained via a force-based stochastic re¬ 
laxation method (described in detail below): All of the 
other nodes (free nodes) are allowed to be displaced un¬ 
til the entire network settles down in a new configuration 
with minimal strain/stress energy minima such that all 
free nodes are in a force balanced state. 


3. FORCE-BASED RELAXATION METHOD 

Here, we employ a stochastic optimization scheme to 
find the force-balanced network configuration under the 
prescribed local perturbations that mimic cell contrac¬ 
tion. The standard stochastic optimization procedure 
consists of a random walk in the configuration space. 
Specifically, in order to find the minimal energy states, 
a randomly selected node of the network is given a ran¬ 
dom small displacement (i.e., a trial move), which leads 
to stretching/compression of fibers that connected to this 
node and thus, causes a change of total elastic energy E G 
of the network. If E G decreases, the displacement is ac¬ 
cepted, otherwise, it is rejected, i.e., the steepest descent 
method is utilized. 

In this scheme, the global energy needs to be re¬ 
computed after each trial displacement. Directly re¬ 
computing E g using Eq. © is both time consuming 
and computationally inefficient due to the large number 
of nodes (5000) and bonds (8000) in our system, as well 
as the large number of trial moves (i.e., ~ 5000 moves 
for each node). Therefore, a local energy update method 
fcrl IyI is employed: before displacing the randomly cho¬ 
sen node, the energy associated with this node, i.e., the 
“local energy” E L , which is defined as the clastic energy 
of all the bonds (fibers) connected to this specific node, 
is calculated. Because the displacement of a node only 
affect the bonds connected to this node, the change of the 
total energy SE G is exactly equal to the change this local 
energy SE L , i.e., SE G = SE L . The acceptance probability 
of the trial move p aC c is simply given by 

r i, se L < o 

Pace = \ (4) 

[ 0, SE^ > 0. 

This also enables of the efficient computation of both the 
global energy and thus, significantly reduces the compu¬ 
tational cost, i.e., from 8 CPU hrs to 0.6 CPU hrs on our 
Dell Precision workstation [with Intel Xeon(R) E5-1603 
2.80GHz 4-core CPU and 16 GB Memory). 

Furthermore, instead of randomly displacing the 
nodes, force-based displacements are used. Specifically, 
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a randomly selected node is displaced along the direction 
of the net force on this node, with the magnitude of the 
displacement proportional to the magnitude of the force. 
For a randomly selected node i, the net force F 1 at this 
node with the components Ff., F* and FI is calculated by 
summing the x , y 1 z components of all the forces exerted 
by all the neighbor nodes on this chosen node i. The 
magnitude of the force between node i and its neighbor 
k connected by a fiber is computed as 

FA 

fik = -jo~- (dik - difc), (5) 

a ik 

where E is the Young’s modulus of the fiber, A is the 
fiber cross section area, dik and d® k are respectively the 
fiber lengths after node displacement and original fiber 
length. Since fik is along the direction defined by nodes 
i and k , the x component of fik can be computed by 
multiplying fik by the x component of the normalized 
direction vector pointing from node i to k, i.e., 


rj-lik 

X 


(■Xk - Xj) 
dik 


( 6 ) 


where Xi and Xk are respectively the x coordinates of 
node i and k. The x component of the net force at node 
i is then given by 


K = E x T *- (?) 

<k> 

The components F* and F\ can be calculated in a similar 
way, which we do not describe in detail here. 

Finally, the node i is displaced along the direction of 
F* with a magnitude 


d° 

S ‘ = C '*' F ‘ < 8 > 

where C is random multiplier between [0, 625,000]. The 
upper bound is chosen such that the maximal individual 
displacement is ~ 1/500 average fiber length, which leads 
to fast convergence of the optimization. About 5000 trial 
moves are performed on each node, with an average suc¬ 
cess rate of 95%. The simulation is terminated when the 
total energy converge to a stable plateau. 

Table I compares the efficiency of different optimiza¬ 
tion schemes described above. Specifically, the different 
procedures yield essentially identical final energy, indi¬ 
cating all of them are robust in finding the desired (lo¬ 
cal) energy minimum. However, the local energy-update 
scheme accelerates the simulation by 12 times, and the 
force-based method possesses a success rate that is 20 
times higher than the random displacement approach. 
The differences in the force-based approach and random- 
displacement approach is also shown in Fig. [3] It can 
be clearly seen that the force-based scheme significantly 
improve convergency of the optimization. 


1. 
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1 

(b) 

FIG. 3: (Color online). Comparison of different optimization 
procedures, (a) The total energy (in pj) as a function of 
trial move number. The red (or dark gray in print version) 
curve corresponds to the force-based method and the black 
curve corresponds to the random displacement method, (b) 
The rate of change of energy. The blue (or dark gray in print 
version) curve corresponds to the force-based method and the 
green (or light gray in print version) curve corresponds to 
the random displacement method. It can be seen that the 
force-based scheme significantly improve convergency of the 
optimization. 




( a ) 


4. SINGLE-CELL CONTRACTION RESULTS 

In this section, the force-based relaxation method is 
applied to study the mechanical response of a biopolymer 
network due to the contraction of a virtual cell embedded 
in the network. In particular, a spherical cell is placed in 
the center of the simulation box, i.e. at (15,15,15 )(pm), 
with a radius R c = (i.e., a linear size of 9 pm). The 

contraction ratio T c (defined as the ratio of the cell radius 
after and before the contraction) is chosen to be 0.3. We 
consider the “cell” is initially residing inside a stress-free 
network. Therefore it is easy for the cell to perturb the 
network via contraction, leading to relatively large T c 
used here. The force-based relaxation method is then 
employed to find force-balanced network configuration. 


4.1. Heterogeneous distribution of forces on fibers 

The distribution of forces on the fibers is shown in 
Fig. [I] It can be clearly seen that the majority of fibers 
(~ 5 000 out of 8 000 fibers in our system) carries very 
small forces (i.e., virtually zero). There are only a small 
number of fibers that carry the majority of the forces 
~ 10 4 pN, which leads to a highly heterogeneous distri¬ 
bution of forces in the network. 


4.2. Identifying force chains 

To understand how the local perturbation due to cell 
contraction is propagated throughout the system, we 
investigate the spatial correlations among large-force¬ 
bearing fibers and identify well-defined structures formed 
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force-based local 

random local 

force-based global 

random global 

tc 

38 min 

40 min 

8 hr 41min 

8 hr 59min 

El 

28.46740 

28.46740 

28.46740 

28.46740 

El 

6.7716672 

6.6983984 

6.7727792 

6.6982896 

N t 

25,000,000 

25,000,000 

25,000,000 

25,000,000 

N s 

23,816,643 

1,143,437 

23,785,917 

1,143,636 

7 

95.3% 

4.57% 

95.1% 

4.57% 


TABLE I: Comparison of different stochastic energy minimization schemes to obtained force-balanced biopolymer network 
configurations. t c is the total CPU computational time consumed, E% is the initial energy (in pj ), E(. is the final energy (in 
pJ), Nt is the total number of MC moves, N s is the number of successfully MC moves, and 7 = N s /Nt is the success rate. 



/ 


FIG. 4: (Color online). Distribution of forces (in pN) on the 
fibers. The majority of fibers (~ 5000 out of 8 000 fiber in 
our system) carries very small forces. There are only a small 
number of fibers that carry the majority of the forces ~ 10 4 
pN. Such fibers form linear structures, which we refer to as 
force chains (see Fig. [5j. 


by these fibers. Specifically, we visualize the entire sys¬ 
tem in a 100 x 100 x 100 voxel-based box, with different 
color code corresponding to different stress states of the 
fibers. The gradual change of color from blue to red indi¬ 
cates the increase of the forces carried by the fibers. Sur¬ 
prisingly, we find that there are well-defined chain-like 
structures composed of large-force-bearing fibers (i.e., 
red fibers) in the network, which are emitting from the 
contracted cell, i.e., the spherical region centered at 
(15,15,15) (/jm) and protruding outwards along the ra¬ 
dial direction, see Fig. [5](a). We refer to these chain-like 
structures as “force chains” in the biopolymer network. 
Figure[5Kb) shows the snapshot of a contracting NIH 3T3 
cell with an almost spherical morphology in collagen I gel 
with a concentration of 1.5mg/ml. Similar linear chain¬ 
like structures of collagen fibers emitting from the cell 
can be clearly observed, which validates our simulation. 

Starting from the contracted spherical cell, we track 
the propagation of forces (generated by the cell) along 
the force chains. In particular, the origin of a force chain 
is associated with a fiber that is directly connected to 
the cell, i.e., with one node as the coarsen-grained “focal 
adhesion” site, points outwardly from the cell and car- 



(b) 


FIG. 5: (Color online). Force chains in the biopolymer net¬ 
work due to the contraction of an embedded spherical cell, 
(a) Simulation result: the fibers carrying large forces are high¬ 
lighted with red color. It can be clearly seen that there are 
chain-like structures emitted from the contracted cell, con¬ 
sisting of fibers bearing very large forces (~ 10 4 pN). (b) 
Experimental validation: A contracting NIH 3T3 cell with an 
almost spherical morphology in collagen I gel with a concen¬ 
tration of 1.5mg/ml. Similar linear chain-like structures of 
collagen fibers emitting from the cell can be clearly observed. 


ries a force larger than a threshold value (e.g., 10 4 pN). 
The next segment of the force chain is then identify as 
the fiber that shares a common node of the origin fiber, 
points outwardly from the cell and carries a force that is 
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larger than all of the other fibers connected to the shared 
node. This process is repeated until the force chain ex¬ 
tends to the boundary of the box. Although it is not the 
only way to identify force chains, our procedure identi¬ 
fies force chains with largest linear extent by imposing 
the condition that the segments of the force chain always 
point outwardly from the cell. This allows us to inves¬ 
tigate the range of fiber-mediated mechanical signaling 
between cells. 


150 

s ioo 

50 

q 

N 

FIG. 6: (Color online). Comparison of the angles a between 
two successive fiber segments along the force chains before 
(black dots) and after (red dots) the cell contraction. This 
comparison reveals two mechanisms for force-chain formation, 
i.e., fiber re-orientation and selection of pre-existing chain-like 
structures in the network, as discussed in the text. 



How does such chain-like structures arise in response to 
the perturbation due to cell contraction? To investigate 
the mechanism for the formation of such linear structure, 
we compare the angles a between two successive fiber 
segments along the force chains, before and after the cell 
contraction. Figure [6] shows the comparison. It can be 
clearly seen that most of the angles a after cell contrac¬ 
tion (shown as red dots) are relatively small (i.e., less 
than 20 degrees), which is consistent with observed lin¬ 
ear chain-like structures. However, a significant fraction 
of a (~ 70%) before contraction (shown as black dots) 
possesses large values (e.g., > 40 degrees). This suggests 
that roughly 70% fibers undergo large re-orientation due 
to cell contraction in order to support the propagation 
of forces through the linear force chains. We note this 
result is consistent with those reported in Ref. 0 - The 
remaining angles, whose values before the contraction are 
relatively small, do not significant change after the per¬ 
turbation. This indicates that certain pre-existing chain¬ 
like structures in the non-stressed network can also be 
selected to carry large forces and contribute to the force 
chains. The relative contributions of the aforementioned 
two mechanisms for force-chain formation, i.e., fiber ori¬ 
entation and selection of pre-existing linear structures, 
generally depends on the network of microstructure. 


4.3. Force decay along force chains 

To quantify how the stress propagates along the force 
chains, the fiber segments along each force chain have 
been identified and the magnitude of the forces on the 
fibers along each force chain are obtained. Figure [7] shows 
the comparison of the decay of forces along the force 
chains fc{r) and the decay of radially averaged force 
f(r) as the distance r form the contracted cell increases. 
We note that f(r) is computed by averaging the magni¬ 
tude of forces carried by the fiber segments whose centers 
are in a concentric thin spherical shell with radius r and 
thickr 000 ' 7 "~ 



FIG. 7: (Color online). Comparison of the decay of forces 
/c(in pN) from the contracting cell along the force chains 
(solid symbols) and the decay of radially averaged force /(in 
pN) in the biopolymer network (dashed line). It can be clearly 
seen that the decay of fc is much slower than decay of /. This 
suggests a possible mechanism for long-range force transmis¬ 
sion in cellularized biopolymer network. 

It can be clearly seen in Fig. [T] that the decay behavior 
of is distinctly different than f(r). In particular, 

/c(r) decays much slower than /(r), although both are 
characterized by a power law decay as can be clearly seen 
in the logarithmic plot. The orange line shows the linear 
fit of In / vs. lnr, i.e., 

ln/(r) = 9.765 - 1.062 Inr. (9) 

The slope 1.062 indicates that the decay of f(r) can be 
very well described via 

f(r) ~ 1 /r, (10) 

which is consistent with the prediction from linear elas¬ 
tic theory (6Cj. Due the large fluctuations of the forces 
along the force chain, the exact decay behavior of fc { r ) 
is difficult to extract, but can be approximated via 

fc(r) ~ 1 /r v , (11) 

where rj e (0.3, 0.5). In addition, the magnitude of the 
forces along the force chain is much higher than the radi¬ 
ally averaged forces. This is consistent with our observa¬ 
tion from the distribution of forces shown in Fig. [2 i.e., 
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the majority of forces are carried by a small number of 
fibers constituting the force chains, while the remaining 
fibers carry very small forces (see fig5 & hg7)0. 


5. DOUBLE-CELL CONTRACTION RESULTS 

To further understand fiber-mediated mechanical cou¬ 
pling between cells, we investigate the force chains 
formed due to the contraction of two cells that are 
close to one another. In particular, two cells with 
R c = 4.5 pm are placed on the diagonal line of the cubic 
simulation box, respectively at (10.5,10.5,10.5 )pm and 
(19.5,19.5,19.5 )pm. We note that these points are the 
trisection points of the body diagonal of the “free” part 
of the simulation box. A contraction ratio of T c = 0.5, 
instead of 0.3 for the single cell case, is employed. This 
is because when more than one cells are embedded in the 
network, the system can be considered in a prestressed 
state and thus, is more difficult to mechanically perturb 
than the single cell case. After the perturbation (i.e., si¬ 
multaneous contraction of both cells), the system is equi¬ 
librated using the force-based relaxation method. Specif¬ 
ically, each node is given on average 7,000 trial moves 
with an averaee success rate of 97.35%. 
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FIG. 8: (Color online). Distribution of the forces (in pN) 
carried by the fibers in the biopolymer network with two con¬ 
tracting cells. It can be seen that the majority of fibers only 
carry very small forces close to zero. The larger forces are 
mainly carried by fibers forming force chaining connecting the 
two contracting cells and emitting outward to the boundary. 


Figure [5] shows the distribution of the forces carried by 
the fibers. As in the single-cell case, the majority of fibers 
only carry very small forces close to zero. The majority 
of forces are carried by a small number of fibers which, 
as shown below, form linear chain-like structures. Figure 
El shows the stressed network due to contraction of the 
cells. The color code is same as in the single-cell case, i.e., 
the fibers carrying large forces are highlighted with red 
color. Two force-chains connecting the two contracting 
cells can be clearly identified, both consisting fiber ele¬ 
ments carrying a forces large than 10 6 pN (see Appendix 



FIG. 9: (Color online). Force chains in the biopolymer net¬ 
work due to the contraction of two embedded cells, (a) Simu¬ 
lation result: The fibers carrying large forces are highlighted 
with red color. It can be clearly seen that there are chain¬ 
like structures connecting the two cells, consisting of fibers 
bearing very large forces (> 10 5 pN). (b) Experimental vali¬ 
dation: Two contracting NIH 3T3 cells close to one another 
in collagen I gel with a concentration of 1.5mg/ml. Similar 
linear chain-like structures of collagen fibers connecting the 
two cells can be clearly observed. 


A). The fibers in the force chains are aligned very well 
with the body diagonal line, on which the two cells are 
placed. Figure Ob) shows the snapshot of two contract¬ 
ing NIH 3T3 cells close to one another in collagen I gel 
with a concentration of 1.5mg/ml. Similar linear chain¬ 
like structures of collagen fibers connecting the two cells 
can be clearly observed, which validates our simulation. 

To understand how such linear-structures form, we 
compute and compare the angles between individual fiber 
segments along the force-chains and the body diagonal, 
before and after cell contraction. Figure [lO] shows the 
comparison. It can be clearly seen that the fiber seg¬ 
ments along the two most prominent force chains between 
the two cells undergo a significant reorientation to align 
with the body diagonal. We note that in this configu¬ 
ration, the largest forces occur along the body diagonal 
connecting the 2 cells, which also drive the re-orientation 
of the fibers. In addition, the fiber segments in the force 










100 

80 

60 

a 

40 
20 

°( 

N 

FIG. 10: (Color online). Comparison the angles between in¬ 
dividual fiber segment along the force-chains and the body 
diagonal, before (black dots) and after (red dots) cell con¬ 
traction. The fiber segments along the force chains undergo 
a significant reorientation to align with the body diagonal, as 
discussed in the text. We note that the angles a shown here 
are with respect to the body diagonal of the simulation box, 
rather than the angles between successive fiber segments in 
the single-cell case. 



chain carrying large forces undergo a more significant re¬ 
orientation (see Appendix A), to get better aligned with 
the body diagonal in order to resist the larger forces. 
The “selection” mechanism for force-chain formation is 
not dominant in the two-cell case. This is because in 
our current network, the probability of finding a chain of 
fibers all perfectly aligned with the body diagonal is very 
small. Finally, we note that since the two-cell configu¬ 
ration does not possess spherical symmetry, we therefore 
do not characterize the radial decay of the forces, which 
is not well defined. The magnitude of forces carried by 
the fiber segments along the force chains are given in 
Appendix A. 


6. CONCLUSIONS AND DISCUSSION 


We have numerically investigated the mechanical re¬ 
sponse of a model biopolymer network due to the ac¬ 
tive contraction of embedded cells. Cell contractions are 
modeled by applying correlated displacements at specific 
nodes, representing the focal adhesion sites. A force- 
based stochastic relaxation method is employed to ob¬ 
tain force-balanced network under cell contraction. We 
find that the majority of the forces are carried by a small 
number of heterogeneous force chains emitted from the 
contracting cells. The force chains consist of fiber seg¬ 
ments that either possess a high degree of alignment be¬ 
fore cell contraction or are aligned due to the reorienta¬ 
tion induced by cell contraction. Large fluctuations of 
the forces along different force chains are observed [c.f. 
Fig. 7]. In addition, the decay of the forces along the 
force chains are significantly slower than the decay of 
concentric averaged forces in the system. These results 


suggest that the fibreous nature of biopolymer network 
structure can support long-range force transmission and 
thus, long-range mechanical signaling between cells [fill ]. 
One of the adavantages of Mechano-signal transduction 
over Chemcial signaling is that it could be upto 40 times 
faster than the diffusion-based, as is reviewed in [fill] . 

We note that force chains are also apparent in gran¬ 
ular materials, which is non-equilibrium state of matter 
composed of macroscopic particles with strong repulsive 
interactions [63j. Compared to the typical force chains in 
granular materials, the force chains we identified in the 
perturbed network have a much better linearity, i.e., they 
either aligned radially outwards in the single-cell case or 
align the body diagonal in the two-cell case. In addition, 
the forces carried by the fiber segments are stretching 
forces, while those carried by the individual grains are 
compressive in nature. The force chains in these two 
distinct systems also share some similar characteristics. 
For example, in both systems, the force chains carry the 
majority of the forces, which are orders of magnitude 
higher than the average forces carried by a fiber/grain. 
This implies that the formation of force chains could be a 
universal mechanism for stress dissipation in disordered 
systems composed a large number of individual building 
blocks with interaction rules in between. 

Also notably, from Figure 7, we can also see the “du- 
rality” of the collagen network system: on one hand, if we 
are interested in the bulk property of the system, it’s to¬ 
tally fine to treat the collagen network as homogeneous in 
the same way as all the other continuum materials, which 
is evident from the perfect 1 /r overall decay of the stress 
field , on the other hand, if we are more intrigued by the 
underlying events happening at even smaller scales, the 
fiberous nature of collagen system can not neglect. 

Finally, we note that our biopolymer network model 
only considers fiber elongation and does not explicitly 
take into account the effect of bending. This assump¬ 
tion works perfectly fine for fibers with short persistent 
lengths, but is not true for fibers with long persistent 
lengths. In future work, the effects of fiber bending will 
be explicitly investigated. In addition, the effects of fiber 
alignment in the original non-stressed network work will 
also be studied. It is expected that in a highly-orientated 
fiber network, the “selection mechanism” will dominant 
for force-chain formations, and thus, significantly biases 
the force propagation as well as the cell migration inside. 
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TABLE II: Realignment of the fibers (i.e., change of the angle 
a (degree) between two successive fiber segments before and 
after contraction) in the force chain between the 2 cells and 
the forces (pN) carried by the fiber segments. 


Element 

a Before 

a After 

Force 


Force Chain 1 


0 

24.87 

2.71 

174,334.2 

1 

21.43 

4.21 

182,937.9 

2 

32.22 

3.83 

182,290.5 

3 

23.85 

3.34 

181,654.4 

4 

38.05 

14.12 

164,516.6 

5 

81.81 

10.51 

152,255.8 

6 

23.30 

12.23 

169,661.7 


Force Chain 2 


0 

20.17 

4.49 

315,013.8 

1 

34.11 

0.81 

219,117.7 

2 

88.69 

0.81 

219,118.4 

3 

23.52 

0.71 

224,975.0 

4 

34.51 

0.72 

224,988.8 

5 

32.98 

0.72 

224,989.5 

6 

61.68 

2.10 

308,953.9 
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